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ABSTRACT 

The precision of radial velocity (RV) measurements depends on the precision attained 
on the wavelength calibration. One of the available options is using atmospheric lines 
as a natural, freely available wavelength reference. Figueira et al. (2010) measured 
the RV of O2 lines using HARPS and showed that the scatter was only of ~10m/s 
over a timescale of 6 yr. Using a simple but physically motivated empirical model, 
they demonstrated a precision of 2m/s, roughly twice the average photon noise con- 
tribution. In this paper we take advantage of a unique opportunity to confirm the 
sensitivity of the telluric absorption lines RV to different atmospheric and observing 
conditions: by means of contemporaneous in-situ wind measurements. This oppor- 
tunity is a result of the work done during site testing and characterization for the 
European Extremely Large Telescope (E-ELT). The HARPS spectrograph was used 
to monitor telluric standards while contemporaneous atmospheric data was collected 
using radiosondes. We quantitatively compare the information recovered by the two 
independent approaches. 

The RV model fitting yielded similar results to that of Figueira et al. (2010), 
with lower wind magnitude values and varied wind direction. The probes confirmed 
the average low wind magnitude and suggested that the average wind direction is a 
function of time as well. However, these results are affected by large uncertainty bars 
that probably result from a complex wind structure as a function of height. The two 
approaches deliver the same results in what concerns wind magnitude and agree on 
wind direction when fitting is done in segments of a couple of hours. Statistical tests 
show that the model provides a good description of the data on all timescales, being 
always preferable to not fitting any atmospheric variation. The smaller the timescale 
on which the fitting can be performed (down to a couple of hours), the better the 
description of the real physical parameters. We conclude then that the two methods 
deliver compatible results, down to better than 5 m/s and less than twice the estimated 
photon noise contribution on O2 lines RV measurement. However, we cannot rule out 
that parameters a and 7 (dependence on airmass and zero-point, respectively) have 
a dependence on time or exhibit some cross-talk with other parameters, an issue 
suggested by some of the results. 

Key words: Atmospheric effects, Instrumentation: spectrographs, Methods: obser- 
vational, Techniques: radial velocities 



1 INTRODUCTION 
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The research on extrasolar planets is currently one of the 
fastest-growing in Astrophysics. Triggered by the pioneering 
work of Mayor & Queloz ( 1995 1 on 51 Peg, it evolved into a 
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domain of its own, with more than 500 planets confirmed up 
to date. Most of these planets (~90%) were detected using 
the radial velocity (RV) induced on the star by the orbital 
motion of the planet around it. The measurement of precise 
RVs can only be done against a precise wavelength reference, 
and two different approaches were pursued extensively. The 
first was the usage of a Th-Ar emission lamp with the cross- 
correlation function (CCF) method (|Baranne et al. 19961, 



and the second the I2 cell along with the deconvolution pro- 
cedure ( Butler et al.|1996 l. In order to measure precise RV 
in the IR with CRIRES ( [Kaeufl et al.|2004| >, |Figueira et al. 
( 2010 1 recovered a method known for a long time: the usage 



of atmospheric features as a wavelength anchor. Using CO2 
lines present in the H band, the authors reached a precision 
of ~5m/s over a timescale of one week. While a similar pre- 
cision had been attained in the past in the optical domain 
using O2 lines, the studies on the stability of atmospheric 
lines were limited to a timescale of up to a couple of weeks. 
In order to assess the RV stability of atmospheric lines 



over longer timescales, Figueira et al. (20101 used HARPS 



(High Accuracy Radial velocity Planet Searcher) archival 
data, spanning more than six years. Three stars - Tau Ceti, 
/i Arae, and e Eri - were selected because they provided a 
strong luminous background against which the atmospheric 
lines could be measured, and were observed not only over a 
long timespan but with high temporal frequency (in astere- 
oseismology campaigns). The spectra were cross-correlated 
against an O2 mask using HARPS pipeline, which delivered 
the RV, bisector span (BIS) and associated uncertainties. 
The high intrinsic stability of HARPS allowed one to mea- 
sure these effects down to lm/s of precision, roughly the 
photon noise attained on the atmospheric lines. 

The r.m.s. of the velocities turned out to be of only 
~10 m/s, and yet well in excess of the attained photon noise. 
An inspection of the RV pattern on one star over one night 
revealed not white noise but a well-defined shape on RV, 
BIS, contrast and FWHM. A component of the RV signal 
was associated with BIS variation, which in turn was lin- 
early correlated with the airmass at which the observation 
was performed. A second component of the signal was in- 
terpreted as being the translation of the atmospheric lines' 
center created by the projection of an average horizontal 
wind vector along the line of sight. These two effects were 
described by the formula 



fi = a x — — - l + /3 x cos{6) . cos{4> - 5) + 7 (1) 
\sm(a) J 

where a is the proportionality constant associated with 
the variation in airmass, f3 and S the average wind speed 
magnitude and direction, and 9 and (f> the telescope elevation 
and azimuth, respectively. The 7 represents the zero-point 
of the RV, which can differ from zero. The fitting of the vari- 
ables a, /3, 7, and 8 allowed a good description of the telluric 
RV signal, with the scatter around the fit being of around 
2 m/s, or twice the photon noise. The fitting was performed 
in two ways: first, allowing all parameters to vary freely and 
second, imposing the same a and 7 for the different datasets. 
For details the interested reader is referred to the original 
paper. However, the model represented by Eq[T] while being 
physically motivated, was not fully validated due to the ab- 
sence of wind measurements against which the fitted values 
could be compared. 



Among the atmospheric parameters studied for E-ELT 
site testing is precipitable water vapor (PWV), the ma- 
jor contributor to the opacity of Earth's atmosphere in 
the infrared. Hence the mean PWV established over long 
timescales determines how well a site is suited for IR as- 
tronomy. For the E-ELT site characterisation a combination 
of remote sensing (satellite data) and on-site data was used 
to derive the mean PWV for several potential sites, taking 



La Silla and Paranal as reference (Kerber et al. 20101. In 



order to better understand the systematics in the archival 
data and to obtain data at higher time resolution, a total of 
three campaigns were conducted at La Silla Paranal obser- 
vatory in 2009. During each campaign all available facility 
instruments as well as dedicated IR radiometers were used 
to measure PWV from the ground (Kerber et al. in prep. 
Querel et al.|20iT l. In addition, radiosondes were launched 
to measure the vertical profile of atmospheric parameters 
in situ, with the goal of calculating the real PWV in the 
atmosphere. Radiosondes are an established standard in at- 
mospheric research and all other methods were validated 
with respect to the radiosonde results with very high fidelity 



( Kerber et al.|2010| [Querel et al.|20To)| Chacon et al.|2010 ) 

In the current paper we present the results of exploiting 
data from the above campaigns: since HARPS observations 
and radiosonde measurements were done in parallel we are 
in a position to make a direct and quantitative compari- 
son of the wind speed parameters (/? and 5). The paper is 
structured as follows. In Sect. 2 we describe the data acquisi- 
tion and reduction of both observing campaigns. Section 3 is 
dedicated to the description of the analysis of data and sub- 
sequent results. In Sect. 4 we discuss the implications of our 
results and we conclude in Sect. 5 with the lessons learned 
from this campaign. 



2 OBSERVATIONS & DATA REDUCTION 



2.1 HARPS measurements 



HARPS (Mayor et al. 20031 is a high-resolution fiber-fed 



cross-dispersed echelle spectrograph installed at the 3.6 m 
telescope at La Silla Observatory. It is characterized by a 
spectral resolution of 110 000 and its 72 orders cover the 
the whole optical range, from 380 to 690 nm. Its extremely 
high stability allows one to measure RV to a precision of bet- 
ter than 60cm/s when a simultaneous Th-Ar lamp is used, 
and of around 1 m/s without the lamp. A dedicated pipeline 
(nicknamed DRS for Data Reduction Software) was created 
to allow for on-the-fly data reduction and RV calculation. 
This pipeline delivers the RV by cross-correlating the ob- 
tained spectra with a weighted binary mask. To calculate 
the atmospheric lines RV variation one needs then only to 
construct a template mask representing the lines to moni- 
tor. This weighted binary mask (|Pepe et al.|2002 1 was built 
using HITRAN database (|Rothman et al.||2005 1 to select 



the O2 lines present in HARPS wavelength domain. For the 
details on HARPS, the data reduction procedure and the 
mask construction, the reader is referred to Figueira ct al. 
(20101. The procedure is identical, with the exception that 



the observations used in the current paper were performed 
without simultaneous Th-Ar. 

For this program, 9 stars were observed: HR3090, 
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HR3476, HR4748, HR5174, HR5987, HR6141, HR6930, 
HR7830, and HR8431, which are fast-rotating A-B stars, 
mostly featureless in the optical domain and suitable to be 
used as telluric standards. For details on the stars the reader 
is referred to t he website "Stars for Meas uring PWV with 
MIKE'Qand to |Thomas-Osip et al.| ( |2007[ ). A total of 1120 
measurements were collected on 8 and 9 of May, 2009, during 
the course of two nights of technical time. The stars were ob- 
served in a complex pattern in such a way that both low and 
high airmass and different patches of the sky were probed 
throughout the night in order to sample any variations of 
PWV. The main consequence is that even a fraction of the 
night with a couple of hours can contain observations of sev- 
eral stars at a wide range of airmass and elevation/azimuth 
coordinates, covering well the independent variables of Eq.[l] 
and allowing a precise estimation of the parameters to be fit. 



2.2 Radiosonde measurements 

The radiosonde (Vaisala RS-92) is a self-contained instru- 
ment package with sensors to measure e.g. temperature and 
humidity combined with a GPS receiver and a radio trans- 
mitter that relays all data in real-time to a receiver on the 
ground. The radiosonde is tied to a helium filled balloon 
and after launch ascends at a rate of a few m/s following 
the prevailing winds. On its ascent trajectory the sonde will 
sample the local atmospheric conditions up to an altitude 
of ~20km, when the balloon will burst. By that time it has 
traveled horizontally ~100 km from the launch site. Since 
it relays its 3D location based on GPS location every two 
seconds, the wind vector exerting force on the balloon can 
be deduced from the change in GPS position. 

A total of 17 radiosondes were launched between the 5 
and the 15 of May of 2009 from La Silla site. One or two 
launches were conducted every day /night. On the 13th no 
data were collected due to a technical problem when radio 
contact with the radiosonde was lost shortly after lunch. 
From the collected physical parameters, the six of interest 
for our study, as well as the nominal precision of the mea- 
surements are presented in Tab. [I] As the sondes rise in 
height, they measure the two horizontal wind components 
on each layer with a nominal precision of lx 10 -3 m/s, much 
higher than that of contemporaneous RV measurements. 

Radiosondes form the backbone of the global net- 
work coordinated by the World Meteorological Organisa- 
tion (WMO) for measuring conditions at the surface and in 
Earth's atmosphere by combining the in-situ atmospheric 
sounding with measurements taken onboard ships aircraft 
and satellites. Coordinated radiosonde launches (one launch 
at 12:00 UTC is the minimum requirement, other launch 
times are 00, 06 and 18 hours UTC) provide a global snap- 
shot of atmospheric conditions which are then used as basic 
input for describing its current state and for modeling future 
conditions. 

The recommended maximum distance between stations 
is 250 km but the global distribution is very uneven and bi- 
ased towards heavily populated areas in the Northern hemi- 
sphere. South America is sparsely covered with Chile oper- 



1 http:/ /www.lco.cl/operations/gmt-site-testing/stars-for- 
measuring-pwv-with-mikc/stars-for-measuring-pwv 



ating 4 stations only one of which (Santo Domingo, WMO 
station number 85586) launches two radiosondes per day 
at 00 and 12 UTC. Data from all active launch sites can be 
found at http://weather.uwyo.edu/upperair/sounding.html 

The WMO also defines the requirements in terms of 
equipment and procedures such as number of barometric 
pressure levels, etc. A number of different radiosondes from 
different manufacturers are used in the various countries. To 
ensure comparability the WMO regularly conducts cross- 
calibration campaigns with parallel measurements (Jauhi- 
ainen & Lehmuskero, 2005Q The Vailsala radiosonde RS- 
92 used in our campaign is considered to be the most reli- 
able and accurate commercial device available. Its minor bi- 
ases in particular for day-time launches are well-documented 
(Jauhiainen & Lehmuskero, 2005, Miloshcvich ct al. 2009). 

The global snapshot of the state of the atmosphere, 
taken at 00:00 and 12:00 UT is used as initial condi- 
tions of global meteorological numerical models (GFS(l), 
ECMWF(2), GME(3) among otherQ. 

These initial conditions are employed in numerical ap- 
proximations using dynamical equations, which predict the 



future state of the atmospheric circulation (Holton 2005) 



The models are a simplification of the atmosphere because 
the horizontal resolution of the grid can be between 60 km 
(GME) to 100 km (GFS) - and sometimes more - and the 
vertical resolution provides only very small number of lay- 
ers, but on a global scale the results are very good and 
have improved considerably over the past decades. There are 
other models called mesoscales models (MM5 (4), WRF(5), 
MesoNH(6), among others), which provide higher spatial 
resolution (horizontal and vertical) use more specific dynam- 
ical equations (physics parameterization) and better resolu- 
tion of the surface terrain. The initial conditions for these 
models are usually the global model augmented by some 
local weather stations. Details on the models are available 
from the sites mentioned above. 

Concerning the applicability of RS data to our purpose 
it is important to note that the radiosonde is the accepted 
standard in atmospheric and meteorological research. For 
global weather forecasting a distance of order 250 km be- 
tween radiosonde launch stations is the desired but by no 
means always achieved standard, while the cadence is be- 



2 Performance of the Vaisala Radiosonde RS92-SGP and 
Vaisala DigiCORA Sounding System MW31 in the WMO 
Mauritius Radiosonde Intercomparison, February 2005. webpage 
http: / / www.vaisala.com/Vaisala %20Documents /White%20Papers / 
Vaisala%20Radiosonde%20RS92%20in%20Mauritius%20Intercomparison.pdf 

3 (1) Global Forecast System jlittp: / / www .emc . ncep .noaa. gov / gmb / 1 
moort hi/gam. html) 

(2) European Center of Medium range Weather Forecasting 

( http : / / www.ecmwf.int /products /data/ operationaLsystem / description / 
brief _history. html) 

(3) Global Numerical Weather Prediction Model 
( |http://journals.ametsoc.org/doi/abs/10.1175/1520-0493(2002) 
130%3C0319%3ATOGIHG%3E2.0.CO%3B2) 

(4) The Fifth-Generation NCAR/ Penn State Mesoscale 
Model. NCAR=National Center of Atmospheric Research 
( jhttp: //www. mmm. near. edu/mm5/| 

(5) Weather Research and Forecasting model (http:/ /www.wrf- 
model.org/index.php ) 

(6) Non-Hydrostatic Mesoscale Atmospheric Model 



I http://mesonh.aero.obs-mip.fr/mesonh/ ) 
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tween 6 and 24 h. Hence, the radiosonde data set that we 
use for comparison with HARPS observations is well within 
the accepted limits of applicability in terms of spatial and 
time resolution. 

It is evident that local topography and diurnal vari- 
ations may limit the value of a set of radiosonde data to 
smaller distances and shorter periods of time. To this end 



there is a very instructive analysis by Kalthoff et al. ( 2002 1 



that is directly applicable to our case. They use the Karl- 
sruhe Atmospheric mesoscale model (KAMM) and compare 
with wind measurements taken at stations around 30 degrees 
South in Chile, including the Cerro Tololo Interamerican 
Observatory (CTIO). La Silla (70°44'4"5 W 29°15T5"4 S) 
is located within that region only about 100 km N of CTIO. 



Kalthoff et al. (2002) find that the wind patterns over this 



region are stable, their diurnal variations are highly repro- 
ducible and that wind conditions are mostly stable during 
night time. Their main finding is that for altitudes between 2 
and 4 km northerly winds prevail whereas above 4 km large 
scale westerly winds dominate. The reason for the Northerly 
wind is a deflection of westerly winds by the Cordillera de 
los Andes which forms a barrier. They provide a physical 
explanation (their section 4) in terms of the Froude number 
(ratio between inertial forces and buoyancy) demonstrating 
that this deflected northerly flow is a naturally stable phe- 
nomenon. As mentioned La Silla is located in the same re- 
gion and the wind roses of Cerro Tololo (2200 m) (their Fig. 
7) and La Silla (2400mJ]are very similar, clearly showing a 
predominance of northerly wind. 

In particular winter months (June-August) are charac- 
terized by very constant daily ground wind properties (Fig 
5 of the same paper). Our observations were made in May. 
In addition it is a well-established fact that wind conditions 
in the free atmosphere are much more stable than in the 



turbulent and highly variable ground layer (see e.g. Holton 
2005l|Wallace fc Hobbs||2006| ). 

As a consequence of the very homogeneous overall wind 
structure between 2 km and 4 km and above 4 km we have 
reason to believe, that the information on the wind vectors 
obtained by a radiosonde will be representative of conditions 
over the time span of at least a good fraction of a night for 
our campaigns. 



3 ANALYSIS & RESULTS 
3.1 HARPS measurements 

We analyzed the data from the 9 stars as if coming from 
a single data set, as there is no reason to treat them sepa- 



rately. As done in Figueira et al. (20101, we discarded the 



27 datapoints with photon noise precision worse than 5 m/s, 
which correspond to only 2.5% of the observations. 

The total RV scatter and average photon noise were 
5.01 m/s and 2.82 m/s, respectively. If one separates the set 
in the two nights that constitute it, the values for the first 
night are of 5.36 and 2.92 m/s, and 4.60 and 2.72 m/s for the 
second night. We note that the photon noise contribution to 
the precision from the stellar spectrum is larger than 1 m/s, 

4 http:/ /www. eso.org/gen-fac/pubs/astclim/lasilla/humidity/ 
LSO_meteo_stat-2002-2006.pdf 



validating the choice of not using the lamp simultaneously 
with the observations. 

We fitted the RV variation on the two nights using Eq.[l] 
as described in Figueira et al. (20101. When fitting, we con- 



sidered splitting the dataset in three different ways and mak- 
ing two different hypothesis for the parameters variation. 
On the splitting of the dataset we employed: 1) the same 
parameters for all the observations, 2) an independent set 
of parameters per night, and 3) a set of parameters for each 
one-third of the night. After allowing all parameters to vary 
freely at first, we repeated this imposing a and 7 to be the 
same for the the whole dataset in 2) and 3). The resulting 
parameters, Xred) an d scatter around the fit are presented 
111 Tab. [2] for each case. The error bars were estimated by 
bootstrapping the residuals and repeating the fitting 10000 
times. The 95% confidence intervals were drawn from the 
distribution of the parameters, and the la uncertainty esti- 
mations are presented. 

While one might be tempted to compare the Xred °f the 
data as a way of quantifying the quality of the fit, there are 
several reasons not to do so. The first is that as one divides 
the data into subsets that are fitted independently, there is 
some ambiguity in how the Xred 01 a se t ls compared with the 
combined Xred °f the subsets. However, more important is 
that we are considering a problem with priors, as the reader 
will realize when noting that /3 € [0, 00 [ . The consequence 
is that this corresponds to the fitting of a non-linear model, 
for which the number of degrees of freedom is ill-defined, 



as recently underlined by Andrae et al. (20101. In order to 



compare the quality of the data description by the different 
models, we follow the recommendations of the same authors. 
We calculate the probability that the normalized residuals 
of the fitting are drawn from a gaussian distribution with 
/j, — and (7 = 1, as expected if no signal is present and the 
scatter is dominated by the measurement uncertainty. To do 
so we use the Kolmogorov-Smirnov test (as implemented in 



Press et al.|1992 1 and compute the probability Pa's which, 
loosely speaking, corresponds to the probability that the 
residuals after fitting the model are drawn from a gaussian 
distribution. The larger the value of Pks, the more appro- 
priate the model is to describe the data-set in hands. We 
also calculated the probability PKs(nofit) for normalized 
RVs of each dataset without fitting the model, but to which 
only the average value was subtracted (which corresponds 
to fitting only a constant). The probability for each case on 
each data set is presented in Tab. [3] 

3.2 Radiosonde measurements 

The measurement of the radiosonde wind vector (it, v) as 
a function of time, or height, while interesting, is hardly 
insightful for our objective. We need to calculate the effect 
of this wind as integrated along the line of sight, such as 
it is measured by any telescope and spectrograph on the 
ground. This will deliver an average wind vector which can 
then be compared with the one obtained with HARPS (see 
the previous section). 

A way of calculating this average wind is to consider 
a plane-parallel atmosphere that is composed of horizontal 
layers. Every radiosonde measurement probes the properties 
of a layer in its ascent. To obtain the average wind speed we 
weight the wind speed of each of these layers with its ab- 
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Parameter Unit Precision [Unit] Comment 



time s 0.1 measurements cadence of one every 2 s 

T K 0.1 

P hPa 0.1 

Height m 1 limited to 30 km 

u m/s 0.01 E-W wind component, East positive 

v m/s 0.01 N-S wind component, North positive 



Table 2. The fitted parameters and data properties, before and after the fitted model is subtracted from it, with all parameters kept 
varying freely (top), and when a and f) are imposed as the same for the datasets (bottom, marked with *). 
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o-(O-C) [m/s] 


cr ph [m/s] 


Y 2 


a [m/s] 


(3 [m/s] 


7 [m/s] 


5[°] 


08+09-05-2009 


1093 


5.01 


4.14 


2.82 


2.20 


7 yo+0.76 
'•' y -0.73 


8 4 7 +0.76 
°'^'-0.68 


220.56io;3i 


126.10tS;£ 


08- 05-2009 (1 st n.) 

09- 05-2009 (2 nd n.) 


554 
539 


5.36 
4.60 


4.38 
3.68 


2.92 
2.72 


2.12 
1.83 


5 74+ 1 - 87 
8-20-0.8? 


1 3 OQ + 312 

fi 75+0-60 
°''°-0.50 


220.32+°;^ 
219.951°;;* 


14 5- 14 ^3 5 00 



(I s * 


u.. 


section 1/3) 


185 


3.46 


2.89 


1.88 


2.25 


8 97+ 1 ' 41 


4 44+ 3 - 31 
^'^-0.89 


223.58lo75 


42 17+ 4 7.29 


(I s * 


u.. 


section 2/3) 


185 


4.58 


4.58 


3.27 


1.71 


10 7Q+ 21 - 27 
±u. I »_x7.gg 


15.99^ 


223.96t49 4 


, 7 99 +77.04 
' -"-96.66 


(I s * 


a., 


section 3/3) 


184 


6.26 


4.51 


3.60 


1.53 


90 77+H.86 
z». I I _ 12 .13 


74.86l«; 9 07 


230.0915:™ 


4.28il 3 3 f 


(2 nd 
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n. 
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08-05-2009 (1 st n.)* 


554 


5.36 


4.38 


2.92 


2.12+ 


6.85 


09-05-2009 (2 nd n.)* 
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4.60 


3.68 


2.72 


1.83+ 


6.85 


global fit parameters 


1093 


5.01 


4.06 


2.82 


1.98 




(I s * n., section 1/3)* 
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3.46 


2.97 


1.88 


2.34+ 


8.81 


(I s * n., section 2/3)* 


185 


4.58 


4.62 


3.27 


1.73+ 
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(1 st n., section 3/3)* 
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6.26 


4.55 
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+1.09 
-1.09 


9.23 


+4.15 
-2.92 


221.39 


+0.42 
-0.39 


- og +36.17 


+1.09 
-1.09 


13.47 


+ 1.82 
-1.59 


221.39 


+0.42 
-0.39 


95.93 


+8.41 
-9.63 



Note that <5 = ° when wind direction points towards North, and positive eastwards. The error bars on each of the fitted parameters 
were drawn by bootstrapping the residuals (see text for details). Note that the X^-ed mar ked with + are not defined in the strict sense: 
they are calculated assuming 4 fitting parameters for the considered subset, with the objective of allowing comparison with the 

corresponding unconstrained fitting. 



sorptivity. In doing so we are considering that the absorption 
line we measure with our spectrograph is the result of the 
product of the transmission of all layers, and each one of 
these creates a small line shifted by its respective horizontal 
wind. It is important to note that we chose doing so because 
absorptivity is proportional to the depth of the line at the 
central wavelength, and thus proportional to the spectral 
information contribution for the CCF as described in |Pepe| 
et all (12002 I. 



The absorptivity on each layer Ai is Ai = 1 - e _T 
where r is the optical depth and calculated as r = I(T) x 
Amplitude l or entz x &o 2 (T, P) x A/i where I is the spectral 



line intensity, Amplitude Lor entz the relative amplitude of a 
Lorentzian function, ao 2 the surface density of O2, and Ah 
the height of the layer in question. 

The first component of r is /, the spectral line in- 
tensity (basically, the line area) and is given in [cm -1 / 
(molecule. cm -2 )] in HITRAN. Since / is a function of T we 
calculated a grid of HITRAN / from the minimum to the 
maximum temperature measured by the radiosondes, with a 
step of 0.1K for all the O2 lines within HARPS wavelength 
domain. For each temperature an average I was assigned to 
the overall spectrum. This gives us I(T), and to obtain val- 
ues for T in between two grid points we fitted second-degree 
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data set 


Pks 


Pjfs(const. fit) 


P KS (no fit) 


08+09-05-2009 


1.47c-12 




4.32e-25 


08-05-2009 (1 st n.) 


2.10e-06 


8.62e-06 


6.28e-19 


09-05-2009 (2 nd n.) 


9.29e-05 


4.44e-04 


2.62e-10 


(1 st n., section 1/3) 


1.36e-01 


4.01c-03 


1.22e-05 


(1 st n., section 2/3) 


4.68e-02 


5.06e-02 


3.02e-02 


(1 st n., section 3/3) 


1.46e-01 


4.22e-02 


1.07e-05 


(2 nd n., section 1/3) 


1.81e-01 


3.31c-02 


5.79e-08 


(2 nd n., section 2/3) 


1.22e-01 


3.41e-01 


4.66e-03 


(2 nd n., section 3/3) 


1.84e-01 


7.93e-02 


5.33e-02 



Table 3. The probability that the residuals and original datasets are drawn from gaussian distributions, as estimated used Kolmogorov- 
Smirnov test (see text for details). We represent by Pfcs t ne datasets in which the parameters can vary freely and Pfcs (const, fit) the 
fit in which a and /3 are imposed to be the same. 



polynomials, which provided a very smooth description of 
the data. Interpolating the values provided the same wind 
values down to 0.01 m/s. 

In order to derive the line depth, one has to apply a cor- 
rection to get the amplitude of the Lorentzian function that 
has the equivalent area, given by 1.0/(n x HWHM). The 
HWHM was set to 1.0, but its absolute value does not affect 
the results significantly, for it affects all layers in the same 
way. Subsequent tests showed that changing it from 0.1 to 10 
led to variations of the order of 0.01 m/s and 0.01 ° on wind 
magnitude and direction, much smaller than the error bars 
of the measurements. The surface density ao 2 (T, P) was cal- 
culated using the ideal gas law and assuming a constant vol- 
ume mixing ratio (VMR) of O2, of 20.946% as function of 
height, which is a reasonable assumption up to 80 km, hence 
well justified in the range of interest of up to 30 km. 

With this we calculated the weighted average of the 
velocities, and a weighted standard deviation as well. The 
error on the average is estimated as being the weighted stan- 
dard deviation. To allow comparison with the results from 
the previous section, one can calculate the vector magni- 
tude and direction. The results are presented in Tab. [4] and 
plotted in Fig. [I] 











1 








' [ 







Figure 1. Evolution of average wind magnitude (upper panel) 
and direction (lower panel), as measured by the radiosondes. The 
two nights during which observations with HARPS were per- 
formed are represented by shadowed, colored zones. The values 
are presented in Tab. [4] 



4 DISCUSSION 
4.1 The RV data 

The first point to note when it comes to the results of Tab. [2] 
is the low value of the r.m.s. of RVs over the two nights: 

5 m/s, less than twice the average photon noise value. As 
one selects smaller sets of data, first individual nights and 
then subsets of these nights, one obtains different fitted pa- 
rameters. This suggests that the parameters are variable on 
a timescale smaller than one night. The higher probability 
Pks presented by the short-timescales datasets, attests to 
the fact that the model is more suitable to test the variabil- 
ity on smaller timescales. The fitting performed imposing 
the same a and 7 provides results with similar quality, bet- 
ter for the complete nights fitting, poorer if the nights are 
divided into subsections. Moreover, when comparing the fit- 
ted data with the raw data, one concludes that fitting the 
data with the models leads to residuals that are closer to 
gaussian than subtracting a constant from it. This means 
that the model description, even when less precise, still de- 



scribes a fraction of the signal contained in the data and is 
always preferable over using raw data. 

It is insightful to compare the data with that of |Figueira] 
|et al.| ( |2010[ ). The fitted a and /3 are comparable, and tend 
to be even lower, while 7 is similar in value. Importantly, S 
varies significantly from one data-set to the other, just like 
it varied between the data sets considered in the previous 
paper. It is important to note that the error bars on some 
measurements are rather large, and the discrepancy between 
some consecutive measurements can be explained by this 
alone. In particular the last two thirds of the first night were 
affected by a high photon noise contribution, and these two 
sets yield the fits with the largest and more asymmetric error 
bars and largest residuals for the unconstrained fit. However, 
the scatter is already smaller than twice the photon noise, 
with or without subtracting the fit. 

While the RV data were obtained with a different scien- 
tific objective than the one presented here (that of determin- 
ing the PWV content in the atmosphere), the observations 
were still done in order to sample as many different patches 
of the sky as possible. However, it is extremely difficult to 
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Table 4. Weighted average of wind components u and v, as well as wind vector magnitude and direction for each of the radiosondes 
launch. 



# probe Observation date and hour u[m/s] v [m/s] ||u + v\\ [m/s] S[°] 



1 


EDT 


/ 


05 


/ 


05 


/ 


09 


/ 


1200 UTC 


-9.43 ± 5.71 


7.91 ± 10.95 


12.31 ± 8.28 


-50.01 ± 42.62 


2 


EDT 


/ 


06 


/ 


05 


/ 


09 


/ 


1200 UTC 


-11.10 ± 5.92 


5.48 ± 8.05 


12.38 ± 6.39 


-63.71 ± 35.54 


Q 
O 


ejlJ 1 


/ 


U I 


/ 


uo 


/ 


no 


/ 


uouu U lO 


o nn -L- 1 c:q 
— y.UU it (.Do 


O Ad -i- 1 
Z.4U it (.ok) 


O Q1 _|_ 1 r;o 
y.ol it l.Oi 


— (O.Ul it 40. UU 


4 


EDT 


/ 


08 


/ 


05 


/ 


09 


/ 


0000 UTC 


3.44 ± 8.55 


4.09 ± 10.27 


5.35 ± 9.60 


40.03 ± 99.69 


5 


EDT 


/ 


08 


/ 


05 


/ 


09 


/ 


0600 UTC 


-0.58 ± 7.76 


7.22 ± 9.54 


7.25 ± 9.53 


-4.58 ± 61.42 


6 


EDT 


/ 


09 


/ 


05 


/ 


09 


/ 


0000 UTC 


-3.43 ± 6.96 


8.61 ± 8.40 


9.27 ± 8.22 


-21.70 ± 44.49 


7 


EDT 


/ 


09 


/ 


05 


/ 


09 


/ 


0600 UTC 


-2.50 ± 7.46 


10.85 ± 8.55 


11.13 ± 8.50 


-12.97 ± 38.69 


8 


EDT 


/ 


09 


/ 


05 


/ 


09 


/ 


1200 UTC 


-3.34 ± 4.82 


16.05 ± 11.78 


16.40 ± 11.57 


-11.75 ± 18.49 


9 


EDT 


/ 


10 


/ 


05 


/ 


09 


/ 


0000 UTC 


5.16 ± 8.98 


9.81 ± 9.48 


11.09 ± 9.37 


27.74 ± 46.99 


10 


EDT 


/ 


10 


/ 


05 


/ 


09 


/ 


0600 UTC 


3.12 ± 6.95 


7.88 ± 7.33 


8.48 ± 7.28 


21.58 ± 47.31 


11 


EDT 


/ 


11 


/ 


05 


/ 


09 


/ 


0000 UTC 


3.24 ± 3.54 


10.67 ± 7.00 


11.15 ± 6.78 


16.89 ± 20.30 


12 


EDT 


/ 


11 


/ 


05 


/ 


09 


/ 


0600 UTC 


2.76 ± 3.34 


10.20 ± 6.31 


10.57 ± 6.16 


15.14 ± 19.63 


13 


EDT 


/ 


12 


/ 


05 


/ 


09 


/ 


0000 UTC 


3.63 ± 3.46 


11.27 ± 6.62 


11.84 ± 6.39 


17.83 ± 18.70 


11 


EDT 


/ 


11 


/ 


05 


/ 


09 


/ 


0000 UTC 


14.69 ± 10.42 


12.54 ± 7.69 


19.31 ± 9.37 


49.53 ± 26.52 


15 


EDT 


/ 


11 


/ 


05 


/ 


09 


/ 


1200 UTC 


11.00 ± 7.69 


8.09 ± 5.96 


13.66 ± 7.13 


53.69 ± 27.77 


16 


EDT 


/ 


15 


/ 


05 


/ 


09 


/ 


0000 UTC 


10.66 ± 8.44 


6.66 ± 4.69 


12.57 ± 7.57 


57.98 ± 27.27 


17 


EDT 


/ 


15 


/ 


05 


/ 


09 


/ 


0600 UTC 


10.32 ± 8.41 


8.84 ± 6.49 


13.59 ± 7.66 


49.43 ± 31.05 



Note that <5 = ° when wind direction points towards North, and positive eastwards. 



4.2 The radiosondes data 

The data from the radiosondes provide some interesting 
clues on the behavior of the atmosphere. The calculation 
of the average horizontal wind was affected by large uncer- 
tainties, a consequence of the complex and variable structure 
of winds as one travels across the atmosphere. However, the 
average wind magnitude is remarkably low, being between 
5 and 16 m/s. The values close in time are in agreement 
within error bars, both in what concerns magnitude and di- 
rection. Note, in particular, the values obtained using probes 
#6,7, and 8, released 6 hours apart and perfectly compatible 
within their assigned uncertainty. It is interesting to point 
out that the measurement yielding the largest uncertainty 
on wind direction is the one with the smallest wind magni- 
tude value, as expected. The wind direction measurements 
suggest also a slow variation of wind direction as a function 
of time. 



4.3 Comparing RV and radiosondes data 

In order to compare the two datasets we computed the time 
center of the observations for each block and selected the 
probe measurement that was closest in time to it. Table [5] 
displays the data in a way that allows an easy comparison of 
the average wind vector magnitude and direction, and the 
difference between the quantities obtained with the two dif- 
ferent methods is presented in Fig.|3] We considered for this 
purpose the unconstrained fit values, for they show higher 
Pks- In what concerns wind magnitude the values from the 
probes agree with the fitted values from RV data, for all 
datasets. The only outlier is the third section of the first 
night, which presents a very large value of (3 and strongly 
asymmetric error bars. As discussed before, the correspond- 
ing HARPS dataset has the largest photon noise contribu- 
tion associated, and very poor azimuth coverage (as can be 
seen in Fig.[2| which can explain the lower quality of the 
fit. In terms of wind direction, the values concerning the 



90 
80 
— 70 

V 60 
3 50 
< 40 
30 

3® 
300 
— 250 

£ 200 
| 150 
< 100 
50 




Altitude and Azimuth as a function of time 





J' 

\ 


v 


V 


7 


s 






V 


\ 




























If 




'J 






p 

















1.5 2.0 
jdb [days] 



2.5 
+ 5.4959e4 



Figure 2. The distribution of elevation and azimuth (8, <j>) as 
function of time, for the two observation nights. The three slices 
used for independent fitting are identified as shadowed regions 
with different colors. 



find suitable stars in all directions and thus sample evenly 
(6, 4>), our independent variables. We succeeded in obtain- 
ing observations for elevation G [30, 85] ° for both nights, 
but most observations were taken between azimuth angles 
4> G [100, 250] °, which might limit the accuracy with which 
the wind direction can be determined. When one looks at 
the distribution of [6, cj>) as a function of time (Fig.|2| one 
concludes that the uneven distribution of the parameters to 
fit can limit the performance of the fit. 

However, the main question that remains is if the model 
can be used imposing the same a and 7 for a given set of 
observations. While the results for the constrained fit are 
slightly worse, they do not allow us to reach a firm conclusion 
with regards to this aspect. 
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| ii 

C 
(0 

■a -100 



dataset [ ] 

Figure 3. Difference between the average wind magnitude (upper 
panel) and direction (lower panel) as measured by the two differ- 
ent methods. The values are presented in Tab. [fj] The full dataset 
fit values are coded in red, those corresponding to single-night 
sets are coded in green, and the night subdivisions are coded in 
blue (electronic version only). 



fit of a subsection of the night agree with those derived 
from radiosonde measurements; those on a longer timescale 
do not. The most straightforward interpretation is that the 
constant horizontal wind hypothesis does not hold for large 
timescales. This is not surprising since wind vectors are vari- 
able over time. In other words, the fit provides a better de- 
scription of the data than no fit, and the residuals are closer 
to gaussianity than the raw data, as seen, but the direction 
has no physical correspondence. However, and as stated be- 
fore, one has to note that the a and a(o-c) values are al- 
ready quite close to the a v u level, and that the ratio between 
any of the former and the latter is smaller than 2. Such ra- 
tio values are smaller than those obtained in Figueira et al. 
( |2010[ ), and we cannot discard the fact that we might be 
approaching the limit of extractable information from the 
current dataset. 

It is arguable that the model might be over-simplistic, 
and the small number of parameters and observables might 
fundamentaly limit the RV signal it can reproduce. Some 
leads point in this direction, and we followed these to pro- 
pose and test alternative models. However, no improvements 
were verified relative to the basic model presented before. We 
present the results of this rather more exploratory digression 
in Appendix A. 

Perfect agreement between the HARPS observations 
and the radiosonde data can not be expected for a number 
of reasons listed below. The HARPS spectrum samples a 
pencil beam through the atmosphere when the star is being 
tracked, while the radiosonde performs in situ-measurements 
along its trajectory governed by the prevailing winds. An- 
other drawback is that the atmosphere is only sampled up 
to an altitude of 20 km; however, at this altitude the density 
of O2 is ten times lower than at the top of the mountain, so 
the weight A% is ten times smaller. Finally, the radiosonde 
is expected to oscillate like a pendulum in its ascent, intro- 
ducing a signal in the measured RV which is not rooted in 



the wind vector it is intended to probe. In spite of all these 
limitations, the two data sets agree and provide a coherent 
picture of the atmospheric impact on RV variation, down to 
better than 5 m/s and less than twice the estimated photon 
noise contribution on O2 lines RV measurement. 

A quantitative assessment of the stability of atmo- 
spheric absorption lines as presented here is of very prac- 
tical value for astronomy. Telluric absorption features are 
imprinted on observations with astronomical spectrographs 
over a wide wavelength range, particularly in the infrared. 
On the one hand this constitutes a complication since the 
features overlay the spectrum of the astronomical target 
leading to blends and line shifts. Hence, all observations 
aiming at a high spectral fidelity in certain regions of the 
infra-red need to correct for atmospheric transmission (e.g. 
|Blake et al.|2010|[Seifahrt et al.|20"To]|Muirhead et al.|2011[ ), 
and not only necessarily in the context of RV measurements 
(e.g 



Uttenthaler et al. 20101. With a comprehensive char- 



acterization of the atmospheric stability, one can assess for 
the first time the impact of considering atmospheric lines to 
be at rest or characterized by a constant speed over a given 
period of time. Their stability (or lack thereof) can explain a 
fraction of the residuals obtained today when fitting the at- 
mospheric transmission with a forward model, which yields 
residuals around a few %. On the other hand the telluric fea- 
tures are also used for wavelength calibration again in partic- 
ular in the infra-red where technical calibration sources are 
less common. For physical reasons atomic spectra emitted by 
lamps show fewer lines and a more uneven distribution than 
in the optical. The Th-Ar hollow cathode lamp is the only 
source whose spectrum has been fully characterized in the 
IR ( |Kerber et al.|2008[ ), and which is being used on CRIRES 
but there are limitations in line density and wavelength cov- 
erage. Gas cells usually also cover only a limited wavelength 
range. As a result telluric absorption lines are an attractive 
alternative in parts of the IR. Based on the results presented 
here the stability of the atmosphere will easily support low 
and medium-reslution spectroscopy, while in particular for 
high resolution and high precision work caution has to be 
applied. The actual wind velocity vector, and its variation, 
during the astronomical observations of course remains un- 
known without independent measurements. Hence, it is not 
possible to derive proper error bars for a quantitative analy- 
sis down to the m/s level unless a full analysis following the 



method described in Figueira et al. (20101 is performed 



5 CONCLUSIONS 

We used HARPS to monitor the RV variation of O2 lines in 
the optical wavelength domain. We compared the fitting of 
a model as described in |Figueira et al.| ( |2010[ ) and the ob- 
tained parameters with those delivered by contemporaneous 
radiosondes. The two approaches deliver the same results in 
what concerns wind magnitude and agree on wind direc- 
tion when fitting is done in chunks of a couple of hours. 
The large uncertainty bars on the values obtained from ra- 
diosondes are likely to be a consequence of a complex wind 
structure as a function of height, a fact that weakens the 
applicability of the assumption of a strong horizontal wind. 
We cannot conclude if the a and 7 parameters should be 
constant as a function of time or not, or if a cross-term be- 
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Table 5. Weighted average of wind components u and v, as well as wind vector magnitude and direction for each of the radiosondes 
launch. 



data set 



P [m/s] 



S[°] 



# probe (time distance) \\u + v\\ [m/s] 



5[°] 



08+09-05-2009 

08- 05-2009 (1 st n.) 

09- 05-2009 (2 nd n.) 
(I s * n., section 1/3) 
(I s * n., section 2/3) 
(I s * n., section 3/3) 
(2 nd n., section 1/3) 
(2 nd n., section 2/3) 
(2 nd n., section 3/3) 



8 47+ u ''° 
°-*'-0.68 

1 3 qq+3 12 



6.75 



+0.60 
-0.50 



4 44+ 3 ' 31 
^■^-0.89 

15.99™ 

74.86+«; ° 

15.98+*;;° 

3 61+ 517 
°-° -0.78 



7.78 



+2.77 
-1.77 



126.10+*$ 

145.14^3 58 0Q 

97 75+ 6 - 67 
at .io_ 7 gl 

42 17+ 47 - 29 

4Z - 1 '-17.60 
i 7 99 +77.04 
1( - zz _96.66 

24 Qfi+ 18 - 08 
90.87^31-05 
66 26+ 17 ' 98 



8 (4h) 
7 (lh) 

10 (1.5k) 

6 (1.5h) 

7 (0.5h) 
7 (2.5k) 

9 (lh) 

10 (lh) 
10 (2k) 



16.40 ± 
11.13 ± 
8.48 ± 
9.27 ± 
11.13 ± 
11.13 ± 
11.09 ± 
8.48 ± 
8.48 ± 



11.57 

8.50 

7.28 

8.22 

8.50 

8.50 

9.37 

7.28 

7.28 



-11.75 ± 
-12.97 ± 
21.58 ± 
-21.70 ± 
-12.97 ± 
-12.97 ± 
27.74 ± 
21.58 ± 
21.58 ± 



18.49 
38.69 
47.31 
44.49 
38.69 
38.69 
46.99 
47.31 
47.31 



Note tkat (5 = 0° when wind direction points towards North, and positive eastwards. 



tween them should be included. However, when these are 
fixed the wind direction does not agree with that extracted 
from the radiosonde, which suggests that the model might be 
incomplete at this level. We tested two different alternative 
models that tried to address this possible incompleteness of 
the physical description, but the results were poorer than 
with the base model. 

Statistical tests showed that the base model provides 
a good description of the data on all timescales, being al- 
ways preferable to not fitting any atmospheric variation, and 
that the smaller the timescale on which it can be performed 
(down to a couple of hours), the better the description of the 
real physical parameters. It is important to note that it is 
for the datasets with higher Pks that the wind parameters 
derived from RV are compatible with those extracted from 
radiosondes measurement. Thus, even though the model pre- 



sented in Figueira et al. ( 2010 1 can probably still be refined, 



the agreement is proven down to better than 5 m/s and less 
than twice the estimated photon noise contribution on O2 
lines RV measurement. 
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APPENDIX A: IMPROVING THE MODEL 

The global picture obtained by analyzing the datasets first 
separately and then together allows one to raise some inter- 



esting questions; within this questions is the potential for 
improving the model. The first point to note is the impact 
on the error bars of the fitted parameters when splitting the 
RV data in subsets for the fitting. When separating the data 
set in two nights, the average error bars for each parameter 
increase by a factor which can be slightly in excess of V2 
(depending on the case, for some the increase is much more 
modest), but when the nights are divided into subsets, the 
increase in the error bars exceeds that expected from the 
reduction of the number of data points used for the fit. In 
particular, the relative increase for the S error bars is much 
larger than those for the other parameters. Another interest- 
ing point comes then into view: the second and third sections 
of the first night show comparatively large error bars for the 
four parameters. However, this point is to be taken carefully 
because, as discussed, the photon noise was higher and the 
azimuth coverage not as complete as for the other datasets. 

These two elements point towards a cross-talk between 
the model parameters. Given the simplicity of the model and 
large number of data points available, it is more likely that 
this behavior comes from trying to fit a too simple model 
rather than being caused by lack of conditioning. 

In addition, a poorer match between wind direction 
from the two methods for the constrained fit suggests that 
a and 7 might not be constant for the datasets at hand. 
This is particularly clear for a, while variations on 7 are 
only of a couple of m/s. Such a behavior can be explained 
by a correlation between a and a parameter which repre- 
sents a quantity expected to change with time. It can also 
be explained if this coefficient has an intrinsic dependence 
on time. And, naturally, it can also be explained by a de- 
pendence of the model parameters - or even the RV itself - 
on a single (unrepresented) quantity. 

The most important hint is probably the high variability 
of 5 and large error bars on its determination: this suggests 
that either the variation associated with this coefficient is 
defined in an incomplete fashion or that some other quantity 
has a similar functional dependence on the parameter which 
the 5 variation tries to accommodate. 

When this information is put together one concludes 
that the most likely improvement to the model is an ad- 
ditional dependence of the airmass impact on RV on the 
direction of the wind. This is not completely unexpected as 
a consequence of the chosen model parametrization. 



InlFigueira et al. ( 2010 \ we had already suggested that if 



the atmosphere has a complex vertical wind structure which 
cannot be represented by a single average wind value, a 
might not be considered as constant. It is so because an 
increased broadening of the CCF (due to the span of veloc- 
ities that displace the absorber) will change the correlation 
coefficient between the broadening and the impact on the 
RV. As a consequence it will change the coefficient between 
airmass and RV, our a. To fully characterize the impact of 
this wind broadening contribution to the a coefficient is ex- 
tremely difficult and requires a line-formation model of the 
atmosphere, which is beyond the scope of this work. How- 
ever, we can propose a refinement of our model in order to 
include this effect, and we test it tentatively in two alterna- 
tive parametrizations to Eq.[T] 



n' = 



sin(p) 



x cos(6) . cos(4> -S) + 7(A1) 
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1 sin{6) ' 



x cos((f>- 5) + t(A2) 



In Eq |Al| we consider a to be dependent on the colin- 
earity with the wind direction. This is expected to be the 
case if there is a scatter of velocity around the central veloc- 
ity p. In this parametrization a contains two components: 
the dependence on airmass and the broadening created by 
the scatter in velocity associated with it. In Eq. |A2| we con- 
sider a variation on this assumption in which only the wind 
direction (and not the projection of this direction along the 
line of sight) has an impact on the measured RV. The fitted 
parameters and quantities associated with each dataset are 
presented in Tab. |A1| The results of the aplication of the 
Kolmogorov-Smirnov test and the Pks derived for the two 
cases are presented in Tab. |A2| 

Unfortunately, these modifications do not lead to an 
improvement. The Pks are smaller than in the previously 
considered cases (and the Xred are larger). The cross-talk 
between the different parameters is increased, with the /3 
parameter reaching zero within the 1-cr uncertainties (given 
the way they were calculated this only means that a large 
number of datasets of the MC was best fitted by /3=0.) As 
a consequence, one is forced to conclude that these alter- 
native models increase the correlation or cross-talk between 
parameters, instead of reducing it. 

One can conceive a model in which the dependence on 
altitude and azimuth is concentrated in the parameters in a 
different way, but this dependence should stem from a phys- 
ical motivation. We are then led to conclude that we prob- 
ably reached the limit of extractable information from this 
dataset and an improvement on the quality of the measure- 
ments is required to take this kind of analysis any further. 
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Table Al. The fitted parameters and data properties, before and after the fitted model is subtracted from it, when Eg, |A1| (top), and 
Eg. |A2| (bottom) are considered. 



data set 


#obs 


a [m/s] 


o-(o-C) [ m / s ] 


<7ph [m/s] 


Y 2 


a [m/s] 


/3 [m/s] 


7 [m/s] 




08+09-05-2009 


1093 


5.01 


4.25 


2.82 


2.35 


22 R7+ 5 ' 34 


7 16+ 1 ' 98 
'■ iD -1.93 


220.5lt°;}» 


15L66l||5 


08- 05-2009 

09- 05-2009 


554 
539 


5.36 
4.60 


4.39 
3.86 


2.92 
2.72 


2.16 
2.09 


1 2 QQ + 7 ' 62 
14 '° 3 -7 .81 

26.96±5;{g 


1541 +4.62 

r> ni +2.83 
°- a -2.38 


219.89t?;o? 
219.94j;g;^ 


156.35^-30 
145. 1313-29 



(I s * 


n.. 


section 1/3) 


185 


3.46 


2.95 


1.88 


2.41 




9 47+5.51 
z -^'-2.47 


999 99+ 1-51 
»■"••"_() .75 


159.31^^49 


(I s * 


n.. 


section 2/3) 


185 


4.58 


4.58 


3.27 


1.71 


-9 49+ 194 11 
».<±y_ 211 . 15 


o o n +97.43 


222 Q4+3- 62 


24.3S™ 3 


(I s * 


n.. 


section 3/3) 


184 


6.26 


4.52 


3.60 


1.54 


25 62+ 33 19 


24 41+34.99 


225.04+4;^ 


27.60l 1 1 5 8 1 9 f 


(2 nd 


ii.. 


, section 1/3) 


180 


4.80 


3.87 


2.52 


2.36 




16 77+ 5 - 29 
1D -' '-16.77 


216.41+. 4 ; 34 


152.53™ 


(2 nd 


ii.. 


, section 2/3) 


180 


4.50 


3.74 


2.84 


1.78 


, , 79 +92.59 
±l.iz_ 75 g7 


2 96+ 5 " 


220.56±J;*? 


71 91 +116.77 
,1 - zi -55.11 


(2 M 


ii.. 


, section 3/3) 


180 


3.86 


3.85 


2.82 


1.91 


11 49+140.58 
_1± -^ z -154.65 


00+ 7 ' 31 




I 79 ^q+58.91 

I I »o»_ 10 7.80 



08+09-05-2009 


1093 


5.01 


4.26 


2.82 


2.37 


20.23l« 2 


4 95+2-37 


220.53 + .g; 2 ° 


151.58i23g 


08-05-2009 


554 


5.36 


4.39 


2.92 


2.18 


10.54l 6 7 ; 98 


14 KQ+5-59 


219. 87+.°;^ 


156.45+. 3 ' 3 ° 


09-05-2009 


539 


4.60 


3.86 


2.72 


2.09 


23.80+J; 73 


1 40+ 3 - 67 


219.97+. [J' 8 , 4 


144.63+J;^ 



(I s * 


n., 


section 1/3) 


185 


3.46 


2.98 


1.88 


2.48 


1 6 1 q+3-51 
1D.10_ 6 g4 


n ro+6.20 


1)91 00+1.27 
zii ' Ji _0 i 84 


159.69t3 5 7 J 3 


(I s * 


n., 


section 2/3) 


185 


4.58 


4.59 


3.27 


1.71 


48 sn+ 148 -22 
-48.80_ 164 41 


16 10+ 10409 


220.82l 7 2 6 73 


153.95tll 3 : f 9 


(I s * 


n., 


section 3/3) 


184 


6.26 


4.52 


3.60 


1.54 


-25.90tg-05 


9Q 8n+ 39 - 14 

za,ou -29 .80 


ooc qq+3.96 
^ z " J -^°— 10.44 


25.391**°° 


(2 nd 


n. 


, section 1/3) 


180 


4.80 


3.91 


2.52 


2.41 


15.60±llf| 


1 R qc+5.78 


216 40+ 4 60 


152.52t 20 5 44 


(2 nd 


n. 


, section 2/3) 


180 


4.50 


3.75 


2.84 


1.78 


Id 94+52.66 
lu - -39 .05 


2 00+ 7 20 
z.uu_ 2 00 


220.43+.1; 74 


91 03+ 84 63 
31 - uo -67 .21 


(2 nd 


n.. 


, section 3/3) 


180 


3.86 


3.83 


2.82 


1.89 


-9 03+ 45 ' 37 
a - uo _54.32 


00+ 6 91 


2l9.0lj:^;fo 


179 46+ 68 ' 26 
i(».w_ 101 71 



Note that <5 = ° when wind direction points towards North, and positive eastwards. The error bars on each of the fitted parameters 

were drawn by bootstrapping the residuals (see text for details). 



data set P KS (n') P KS (tt") P K s(no fit) 



08+09-05-2009 


2.51c- 


-11 


6. 


,74c- 


-11 


4.32e- 


-25 


08-05-2009 


7.64b- 


-07 


3 


19e- 


-07 


6.28b- 


■19 


09-05-2009 


3.08e- 


■06 


3 


.51e- 


-06 


2.62b- 


■10 


(I s * n., section 1/3) 


1.15b- 


•02 


8. 


,14e- 


•03 


1.22b- 


-05 


(I s * n., section2/3) 


4.17b- 


•02 


3 


.22e- 


■02 


3.02b- 


-02 


(1 st n., section 3/3) 


1.66b- 


■01 


1 


,84e- 


-01 


1.07e- 


-05 


(2 nd n., section 1/3) 


1.09b- 


•03 


1 


36e- 


■03 


5.79b- 


-08 


(2 nd n., section 2/3) 


1.03b- 


■01 


9. 


,67e- 


■02 


4.66e- 


•03 


(2 nd n., section 3/3) 


2.95b- 


■02 


4. 


,80e- 


■02 


5.33b- 


-02 



Table A2. The probability that the residuals and original datasets are drawn from gaussian distributions, as estimated used Kolmogorov- 
Smirnov test (see text for details), now for the alternative models represented by Eg. |A1| (PjysffiQ), and Eg. |A2| (Pyg(ti")). 



